################################################################################
### "Someone like me? Disability identity and representation perceptions",
### Stefanie Reher and Elizabeth Evans, Political Behavior
### 24 July 2024
################################################################################

################################################################################
### Preparations ###############################################################
################################################################################

rm(list=ls())
options(scipen=999)

### Install and load necessary packages 

# If a package is not yet installed, remove hash to install it

#install.packages("clubSandwich")
library(clubSandwich)
#install.packages("lme4")
library(lme4)
#install.packages("lmerTest")
library(lmerTest)
#install.packages("plm")
library(plm)
#install.packages("lmtest")
library(lmtest)
#install.packages("ggplot2")
library(ggplot2)
#install.packages("ggpubr")
library(ggpubr)
#install.packages("dplyr")
library(dplyr)
#install.packages("mediation")
library(mediation)
#install.packages("sjPlot")
library(sjPlot)
#install.packages("sjmisc")
library(sjmisc)
#install.packages("effects")
library(effects)
#install.packages("stargazer")
library(stargazer)
#install.packages("gmodels")
library(gmodels)


### Create function to generate data frame with estimates from coeftest to create coefficient plots

ctdf <- function(x){
  rt=list()                             # generate empty results list
  for(c in 1:dim(x)[2]) rt[[c]]=x[,c]   # writes column values of x to list
  rt=as.data.frame(rt)                  # converts list to data frame object
  names(rt)=names(x[1,])                # assign correct column names
  rt[,"sig"]=symnum(rt$`Pr(>|z|)`, corr = FALSE, na = FALSE,
                    cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),symbols = c("***", "**", "*", ".", " "))
  return(rt)
}

### Set working directory and load data

# setwd: add your own working directory
setwd(" ")
data <- read.csv("data.csv")

### Assign reference categories of factor variables

data$disid3 <- relevel(as.factor(data$disid3), ref = "Non-disabled")
data$c.dis <- relevel(as.factor(data$c.dis), ref="None")

### Set seed for mediation analysis:
set.seed(1000)

################################################################################
### Analysis ###################################################################
################################################################################

##### Construct preference congruence measure: mean(1 - abs(candidate preference - respondent preference))

data$cong1 <- 1 - abs(data$imp1.nor - data$imp.own1.n)
data$cong2 <- 1 - abs(data$imp2.nor - data$imp.own2.n)
data$cong3 <- 1 - abs(data$imp3.nor - data$imp.own3.n)
data$cong4 <- 1 - abs(data$imp4.nor - data$imp.own4.n)
data$cong5 <- 1 - abs(data$imp5.nor - data$imp.own5.n)
data$cong6 <- 1 - abs(data$imp6.nor - data$imp.own6.n)
data$conglr <- 1 - abs(data$lr.s - data$lrown.s)

data <- data %>%  rowwise() %>%
  mutate(congboth = mean(c_across(c('cong1', 'cong2', 'cong3', 'cong4', 'cong5', 'cong6', 'conglr')), na.rm=TRUE))
summary(data$congboth) 

################################################################################
################################################################################
################################################################################

##### Table 1. Frequencies of disability, disability types, and disability identity in the samples

### UK:

# Disability: 
table(data$disown [data$country=="UK"], useNA = "a")
# total number of disabled respondents: 2952/4 = 738, since there are 4 observations per respondent
2952 / (2952 + 8728) # 25.3% disabled (NA excluded from total)

# Disability identity:
table(data$disid [data$country=="UK"], useNA = "a")
# total number of disability identifiers: 756/4 = 189
756 / (1876 + 756) # 28.7% of disabled people have disability ID
756 / (1876 + 756 + 8728) # 6.7% of all respondents have disability ID

### US:

# Disability:
table(data$disown [data$country=="US"], useNA = "a") 
# total number of disabled respondents: 3448/4 = 862, since there are 4 observations per respondent
3448 / (3448 + 8304) # 29.3% disabled 

# Disability identity:
table(data$disid [data$country=="US"], useNA = "a")
# total number of disability identifiers: 1192/4 = 298
1192 / (1192 + 1844) # 39.3% of disabled people have disability ID
1192 / (1192 + 1844 + 8304) # 10.5% of all respondents have disability ID

################################################################################
################################################################################
################################################################################

#### Supplementary Information 2, Table S2: Linear regression of representation 
# perceptions on candidate disability 

### Model 1: no interaction with disability identity

rep1.m1 <- plm(rep1.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", 
               data=subset(data, !is.na(disid3)))
rep1.m1.coef <- coeftest(rep1.m1, vcov=vcovHC(rep1.m1, type="HC0", cluster="group"))  
rep1.m1.coef

################################################################################

### Model 2 (shown in Figures 2 and 3)

rep1.m2 <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="Non-disabled") + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", 
               data)
rep1.m2.coef <- coeftest(rep1.m2, vcov=vcovHC(rep1.m2, type="HC0", cluster="group"))  
rep1.m2.coef
rep1.m2.frame=ctdf(rep1.m2.coef)

# same model but changed reference categories for disability identity (for Figures 2 and 3)

rep1.m3 <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="No ID") + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", 
               data)
rep1.m3.coef <- coeftest(rep1.m3, vcov=vcovHC(rep1.m3, type="HC0", cluster="group"))  
rep1.m3.coef
rep1.m3.frame=ctdf(rep1.m3.coef)

rep1.m4 <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="ID") + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", 
               data)
rep1.m4.coef <- coeftest(rep1.m4, vcov=vcovHC(rep1.m4, type="HC0", cluster="group"))  
rep1.m4.coef
rep1.m4.frame=ctdf(rep1.m4.coef)

################################################################################

### Figure 2: Effects of candidate disability on representation perceptions, by disability identity 

labels <- c("Respondent", "Coefficient", "SE")
row1 <- c("Non-disabled", rep1.m2.frame[2,1], rep1.m2.frame[2,2])
row2 <- c("No ID", rep1.m3.frame[2,1], rep1.m3.frame[2,2])
row3 <- c("ID", rep1.m4.frame[2,1], rep1.m4.frame[2,2])

all <- data.frame(rbind(row1,row2,row3))
names(all) <- labels
interval <- 1.96 
all$Coefficient <- as.numeric(as.character(all$Coefficient))
all$SE <- as.numeric(as.character(all$SE))
all$Respondent <- factor(all$Respondent, levels=c("Non-disabled", "No ID", "ID"))

plot <- ggplot(all)
plot <- plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
plot <- plot + geom_pointrange(aes(x = Respondent, y = Coefficient, ymin = Coefficient - SE*interval,
                                   ymax = Coefficient + SE*interval),   fill = "WHITE")
plot <- plot +  theme(plot.title = element_text(hjust = 0.5), 
                      axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0))) + 
  theme_bw() + ggtitle(" ") + xlab("Voter disability identity") + coord_flip() + ylim(-0.1, 0.1) +
  ylab("Effect of candidate disability on representation perceptions") 
print(plot)
ggsave("Fig2.png", dpi=300)

################################################################################

### Figure 3: Predicted representation perceptions of disabled and non-disabled 
# candidates, by disability identity 

data$candidate.f <- as.factor(data$candidate) # can't have as.factor within the glm() function to use ggeffects
data$Disability <- NA
data$Disability[data$c.disdummy==0] <- "Not disabled"
data$Disability[data$c.disdummy==1] <- "Disabled"

data$disid3 <- ordered(data$disid3, levels=c("Non-disabled", "No ID", "ID"))

rep1.m2.pred <- glm(rep1.nor ~ Disability*disid3 + c.age + c.female + c.minority + 
                      c.job + c.exp + c.child.num + c.heldoffice + c.party3 + 
                      country + candidate.f, 
                    data=data)

plot <- plot_model(rep1.m2.pred, type = "eff", terms = c("disid3", "Disability"), 
                   dot.size = 2, ci.style = c("bar"), dodge = 0.5) +
  xlab("Voter disability identity") +  ylab("Predicted representation perception") + 
  ggtitle(" ")  +
  theme_bw() + ylim(0.5, 0.7) + coord_flip() +
  scale_colour_manual(values = c("black", "gray60")) + 
  labs(color = "Candidate disability")
plot
ggsave("Fig3.png", dpi=300)

################################################################################
################################################################################

### Table S2, Model 3 (UK only)

# recode disid3 into unordered factor
data$disid3 <- as.character(data$disid3)
data$disid3 <- as.factor(data$disid3)

rep1.m2.uk <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="Non-disabled") + c.age + 
                    c.female + c.minority + c.job + c.exp + c.child.num + 
                    c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                  data=subset(data, country=="UK"))
rep1.m2.uk.coef <- coeftest(rep1.m2.uk, vcov=vcovHC(rep1.m2.uk, type="HC0", cluster="group"))  
rep1.m2.uk.coef
rep1.m2.uk.frame=ctdf(rep1.m2.uk.coef)

# same model with other reference categories of disability identity for Figure S2

rep1.m3.uk <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="No ID") + c.age + 
                    c.female + c.minority + c.job + c.exp + c.child.num + 
                    c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                  data=subset(data, country=="UK"))
rep1.m3.uk.coef <- coeftest(rep1.m3.uk, vcov=vcovHC(rep1.m3.uk, type="HC0", cluster="group"))  
rep1.m3.uk.coef
rep1.m3.uk.frame=ctdf(rep1.m3.uk.coef)

rep1.m4.uk <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="ID") + c.age + 
                    c.female + c.minority + c.job + c.exp + c.child.num + 
                    c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                  data=subset(data, country=="UK"))
rep1.m4.uk.coef <- coeftest(rep1.m4.uk, vcov=vcovHC(rep1.m4.uk, type="HC0", cluster="group"))  
rep1.m4.uk.coef
rep1.m4.uk.frame=ctdf(rep1.m4.uk.coef)

################################################################################

### Table S2, Model 4 (US only)

rep1.m2.us <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="Non-disabled") + c.age + 
                    c.female + c.minority + c.job + c.exp + c.child.num + 
                    c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                  data=subset(data, country=="US"))
rep1.m2.us.coef <- coeftest(rep1.m2.us, vcov=vcovHC(rep1.m2.us, type="HC0", cluster="group"))  
rep1.m2.us.coef
rep1.m2.us.frame=ctdf(rep1.m2.us.coef)

# same model with other reference categories of disability identity for Figure S2

rep1.m3.us <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="No ID") + c.age + 
                    c.female + c.minority + c.job + c.exp + c.child.num + 
                    c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                  data=subset(data, country=="US"))
rep1.m3.us.coef <- coeftest(rep1.m3.us, vcov=vcovHC(rep1.m3.us, type="HC0", cluster="group"))  
rep1.m3.us.coef
rep1.m3.us.frame=ctdf(rep1.m3.us.coef)

rep1.m4.us <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="ID") + c.age + 
                    c.female + c.minority + c.job + c.exp + c.child.num + 
                    c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                  data=subset(data, country=="US"))
rep1.m4.us.coef <- coeftest(rep1.m4.us, vcov=vcovHC(rep1.m4.us, type="HC0", cluster="group"))  
rep1.m4.us.coef
rep1.m4.us.frame=ctdf(rep1.m4.us.coef)

################################################################################

### Figure S2: 

# UK:

labels <- c("Respondent", "Coefficient", "SE")
row1 <- c("Non-disabled", rep1.m2.uk.frame[2,1], rep1.m2.uk.frame[2,2])
row2 <- c("No ID", rep1.m3.uk.frame[2,1], rep1.m3.uk.frame[2,2])
row3 <- c("ID", rep1.m4.uk.frame[2,1], rep1.m4.uk.frame[2,2])

all <- data.frame(rbind(row1,row2,row3))
names(all) <- labels
all$Coefficient <- as.numeric(as.character(all$Coefficient))
all$SE <- as.numeric(as.character(all$SE))
all$Respondent <- factor(all$Respondent, levels=c("Non-disabled", "No ID", "ID"))

plot <- ggplot(all)
plot <- plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
plot <- plot + geom_pointrange(aes(x = Respondent, y = Coefficient, ymin = Coefficient - SE*interval,
                                   ymax = Coefficient + SE*interval),   fill = "WHITE")
plot <- plot +  theme(plot.title = element_text(hjust = 0.5), 
                      axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0))) + 
  theme_bw() + ggtitle("                                   UK") + xlab("Voter disability identity") + coord_flip() + ylim(-0.1, 0.103) +
  ylab("Effect of candidate disability on representation perceptions") 
print(plot)
ggsave("FigS2-UK.png", dpi=300)

# US:

labels <- c("Respondent", "Coefficient", "SE")
row1 <- c("Non-disabled", rep1.m2.us.frame[2,1], rep1.m2.us.frame[2,2])
row2 <- c("No ID", rep1.m3.us.frame[2,1], rep1.m3.us.frame[2,2])
row3 <- c("ID", rep1.m4.us.frame[2,1], rep1.m4.us.frame[2,2])

all <- data.frame(rbind(row1,row2,row3))
names(all) <- labels
all$Coefficient <- as.numeric(as.character(all$Coefficient))
all$SE <- as.numeric(as.character(all$SE))
all$Respondent <- factor(all$Respondent, levels=c("Non-disabled", "No ID", "ID"))

plot <- ggplot(all)
plot <- plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
plot <- plot + geom_pointrange(aes(x = Respondent, y = Coefficient, ymin = Coefficient - SE*interval,
                                   ymax = Coefficient + SE*interval),   fill = "WHITE")
plot <- plot +  theme(plot.title = element_text(hjust = 0.5), 
                      axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0))) + 
  theme_bw() + ggtitle("                                   US") + xlab("Voter disability identity") + coord_flip() + ylim(-0.1, 0.103) +
  ylab("Effect of candidate disability on representation perceptions") 
print(plot)
ggsave("FigS2-US.png", dpi=300)

################################################################################
################################################################################

### Table S2, Model 5 (by disability type)

rep1.m2.dis <- plm(rep1.nor ~ c.dis*relevel(disid3, ref="Non-disabled") + c.age + c.female + c.minority + c.job + c.exp +
                     c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID",
                   data)
rep1.m2.dis.coef <- coeftest(rep1.m2.dis, vcov=vcovHC(rep1.m2.dis, type="HC0", cluster="group"))
rep1.m2.dis.coef
rep1.m2.dis.frame=ctdf(rep1.m2.dis.coef)

rep1.m3.dis <- plm(rep1.nor ~ c.dis*relevel(disid3, ref="No ID") + c.age + c.female + c.minority + c.job + c.exp +
                     c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID",
                   data)
rep1.m3.dis.coef <- coeftest(rep1.m3.dis, vcov=vcovHC(rep1.m3.dis, type="HC0", cluster="group"))
rep1.m3.dis.coef
rep1.m3.dis.frame=ctdf(rep1.m3.dis.coef)

rep1.m4.dis <- plm(rep1.nor ~ c.dis*relevel(disid3, ref="ID") + c.age + c.female + c.minority + c.job + c.exp +
                     c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID",
                   data)
rep1.m4.dis.coef <- coeftest(rep1.m4.dis, vcov=vcovHC(rep1.m4.dis, type="HC0", cluster="group"))
rep1.m4.dis.coef
rep1.m4.dis.frame=ctdf(rep1.m4.dis.coef)

################################################################################

### Figure S3: Effects of candidate disability types on representation perceptions, 
# by disability identity 

labels <- c("Disability", "Respondent", "Coefficient", "SE")
row1 <- c("Blind", "Non-disabled", rep1.m2.dis.frame[2,1], rep1.m2.dis.frame[2,2])
row2 <- c("Deaf", "Non-disabled", rep1.m2.dis.frame[3,1], rep1.m2.dis.frame[3,2])
row3 <- c("Wheelchair user", "Non-disabled", rep1.m2.dis.frame[4,1], rep1.m2.dis.frame[4,2])
row4 <- c("Blind", "No ID", rep1.m3.dis.frame[2,1], rep1.m3.dis.frame[2,2])
row5 <- c("Deaf", "No ID", rep1.m3.dis.frame[3,1], rep1.m3.dis.frame[3,2])
row6 <- c("Wheelchair user", "No ID", rep1.m3.dis.frame[4,1], rep1.m3.dis.frame[4,2])
row7 <- c("Blind", "ID", rep1.m4.dis.frame[2,1], rep1.m4.dis.frame[2,2])
row8 <- c("Deaf", "ID", rep1.m4.dis.frame[3,1], rep1.m4.dis.frame[3,2])
row9 <- c("Wheelchair user", "ID", rep1.m4.dis.frame[4,1], rep1.m4.dis.frame[4,2])

all <- data.frame(rbind(row1,row2,row3,row4,row5,row6,row7,row8,row9))
names(all) <- labels
all$Coefficient <- as.numeric(as.character(all$Coefficient))
all$SE <- as.numeric(as.character(all$SE))
all$Respondent <- factor(all$Respondent, levels=c("Non-disabled", "No ID", "ID"))
all$Disability <- factor(all$Disability, levels=c("Wheelchair user","Deaf","Blind" ))

plot <- ggplot(all)
plot <- plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
plot <- plot + geom_pointrange(aes(x = Respondent, y = Coefficient, ymin = Coefficient - SE*interval,
                                   ymax = Coefficient + SE*interval, colour=Disability, shape=Disability),
                               lwd = 1/2, position = position_dodge(width = 1/2),
                               fill = "WHITE")
plot <- plot +  theme(plot.title = element_text(hjust = 0.5),
                      axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0))) +
  theme_bw() + ggtitle(" ") + xlab("Voter disability identity") + coord_flip() + ylim(-0.1, 0.1)
print(plot + scale_colour_manual(name="Candidate disability \n(ref=Non-disabled)", 
                                 values = c("black", "black", "black")) +
        ylab("Effect of candidate disability on representation perceptions") +   
        scale_shape_manual(guide="none", values=c(15,16,17)) +
        guides(colour = guide_legend(override.aes = list(shape = c(17,16,15)), reverse = T)))
ggsave("FigS3.png", dpi=300)


################################################################################
################################################################################

### Table S2 

stargazer(rep1.m1.coef, rep1.m2.coef, rep1.m2.uk.coef, 
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS2a.doc")

stargazer(rep1.m2.us.coef, rep1.m2.dis.coef,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS2b.doc")

# to get the N for each model: 
stargazer(rep1.m1, rep1.m2, rep1.m2.uk, rep1.m2.us, rep1.m2.dis,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS2-n.doc")

################################################################################
################################################################################
################################################################################

### Figure S4: Effects of candidate disability types on representation perceptions, 
# by disability identity, excluding respondents with and without mobility impairments

# Variable disown.1: 1 if respondents ticked category "Mobility", 0 if not

data$disid3 <- as.character(data$disid3)
data$disid3 <- as.factor(data$disid3)

# Proportion of respondents with mobility impairment
table(data$disown.1, data$disown)
2480 / (2480 + 3920) # 38.8% of disabled respondents have a mobility impairment

### Only including disabled respondents *without* mobility impairments

rep1.m2.nmob <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="Non-disabled") + c.age + c.female + c.minority + c.job + c.exp + 
                      c.child.num + c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                    data=subset(data, disown.1!=1))
rep1.m2.nmob.coef <- coeftest(rep1.m2.nmob, vcov=vcovHC(rep1.m2.nmob, type="HC0", cluster="group"))  
rep1.m2.nmob.coef
rep1.m2.nmob.frame=ctdf(rep1.m2.nmob.coef)

rep1.m3.nmob <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="No ID") + c.age + c.female + c.minority + c.job + c.exp + 
                      c.child.num + c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                    data=subset(data, disown.1!=1))
rep1.m3.nmob.coef <- coeftest(rep1.m3.nmob, vcov=vcovHC(rep1.m3.nmob, type="HC0", cluster="group"))  
rep1.m3.nmob.coef
rep1.m3.nmob.frame=ctdf(rep1.m3.nmob.coef)

rep1.m4.nmob <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="ID") + c.age + c.female + c.minority + c.job + c.exp + 
                      c.child.num + c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                    data=subset(data, disown.1!=1))
rep1.m4.nmob.coef <- coeftest(rep1.m4.nmob, vcov=vcovHC(rep1.m4.nmob, type="HC0", cluster="group"))  
rep1.m4.nmob.coef
rep1.m4.nmob.frame=ctdf(rep1.m4.nmob.coef)

### Only including disabled respondents *with* mobility impairments

rep1.m2.mob <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="Non-disabled") + c.age + c.female + c.minority + c.job + c.exp + 
                     c.child.num + c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                   data=subset(data, (disown.1==1) | (disown=="Non-disabled")))
rep1.m2.mob.coef <- coeftest(rep1.m2.mob, vcov=vcovHC(rep1.m2.mob, type="HC0", cluster="group"))  
rep1.m2.mob.coef
rep1.m2.mob.frame=ctdf(rep1.m2.mob.coef)

rep1.m3.mob <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="No ID") + c.age + c.female + c.minority + c.job + c.exp + 
                     c.child.num + c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                   data=subset(data, (disown.1==1) | (disown=="Non-disabled")))
rep1.m3.mob.coef <- coeftest(rep1.m3.mob, vcov=vcovHC(rep1.m3.mob, type="HC0", cluster="group"))  
rep1.m3.mob.coef
rep1.m3.mob.frame=ctdf(rep1.m3.mob.coef)

rep1.m4.mob <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="ID") + c.age + c.female + c.minority + c.job + c.exp + 
                     c.child.num + c.heldoffice + c.party3 + as.factor(candidate), model="p", index="ID", 
                   data=subset(data, (disown.1==1) | (disown=="Non-disabled")))
rep1.m4.mob.coef <- coeftest(rep1.m4.mob, vcov=vcovHC(rep1.m4.mob, type="HC0", cluster="group"))  
rep1.m4.mob.coef
rep1.m4.mob.frame=ctdf(rep1.m4.mob.coef)

################################################################################

## Figure S4:

# Left figure: only not mobility impaired

labels <- c("Respondent", "Coefficient", "SE")
row1 <- c("Non-disabled", rep1.m2.nmob.frame[2,1], rep1.m2.nmob.frame[2,2])
row2 <- c("No ID", rep1.m3.nmob.frame[2,1], rep1.m3.nmob.frame[2,2])
row3 <- c("ID", rep1.m4.nmob.frame[2,1], rep1.m4.nmob.frame[2,2])

all <- data.frame(rbind(row1,row2,row3))
names(all) <- labels
interval <- 1.96 
all$Coefficient <- as.numeric(as.character(all$Coefficient))
all$SE <- as.numeric(as.character(all$SE))
all$Respondent <- factor(all$Respondent, levels=c("Non-disabled", "No ID", "ID"))

plot <- ggplot(all)
plot <- plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
plot <- plot + geom_pointrange(aes(x = Respondent, y = Coefficient, ymin = Coefficient - SE*interval,
                                   ymax = Coefficient + SE*interval),   fill = "WHITE")
plot <- plot +  theme(plot.title = element_text(hjust = 0.5), 
                      axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0))) + 
  theme_bw() + ggtitle("Not mobility impaired") + xlab("Voter disability identity") + coord_flip() + ylim(-0.1, 0.103) +
  ylab("Effect of candidate disability on representation perceptions") 
print(plot)
ggsave("FigS4-left.png", dpi=300)

# Right figure: only mobility impaired

labels <- c("Respondent", "Coefficient", "SE")
row1 <- c("Non-disabled", rep1.m2.mob.frame[2,1], rep1.m2.mob.frame[2,2])
row2 <- c("No ID", rep1.m3.mob.frame[2,1], rep1.m3.mob.frame[2,2])
row3 <- c("ID", rep1.m4.mob.frame[2,1], rep1.m4.mob.frame[2,2])

all <- data.frame(rbind(row1,row2,row3))
names(all) <- labels
all$Coefficient <- as.numeric(as.character(all$Coefficient))
all$SE <- as.numeric(as.character(all$SE))
all$Respondent <- factor(all$Respondent, levels=c("Non-disabled", "No ID", "ID"))

plot <- ggplot(all)
plot <- plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
plot <- plot + geom_pointrange(aes(x = Respondent, y = Coefficient, ymin = Coefficient - SE*interval,
                                   ymax = Coefficient + SE*interval),   fill = "WHITE")
plot <- plot +  theme(plot.title = element_text(hjust = 0.5), 
                      axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0))) + 
  theme_bw() + ggtitle("Only mobility impaired") + xlab("Voter disability identity") + coord_flip() + ylim(-0.1, 0.103) +
  ylab("Effect of candidate disability on representation perceptions") 
print(plot)
ggsave("FigS4-right.png", dpi=300)

################################################################################
################################################################################
################################################################################

#### SI 3: Disability effects on perceived candidate preferences and citizen preferences

###  Table S3 & Figure S5: Linear regression of perceived candidate preferences on candidate disability

### Social security/welfare 
imp1.m1 <- plm(imp1.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", data)
imp1.m1.coef <- coeftest(imp1.m1, vcov=vcovHC(imp1.m1, type="HC0", cluster="group"))  
imp1.m1.coef
imp1.m1.frame=ctdf(imp1.m1.coef) 

### Military & defence
imp2.m1 <- plm(imp2.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", data)
imp2.m1.coef <- coeftest(imp2.m1, vcov=vcovHC(imp2.m1, type="HC0", cluster="group"))  
imp2.m1.coef
imp2.m1.frame=ctdf(imp2.m1.coef)

### Healthcare
imp3.m1 <- plm(imp3.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", data)
imp3.m1.coef <- coeftest(imp3.m1, vcov=vcovHC(imp3.m1, type="HC0", cluster="group"))  
imp3.m1.coef
imp3.m1.frame=ctdf(imp3.m1.coef)

### Minority rights
imp4.m1 <- plm(imp4.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", data)
imp4.m1.coef <- coeftest(imp4.m1, vcov=vcovHC(imp4.m1, type="HC0", cluster="group"))  
imp4.m1.coef
imp4.m1.frame=ctdf(imp4.m1.coef)

### Economy
imp5.m1 <- plm(imp5.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", data)
imp5.m1.coef <- coeftest(imp5.m1, vcov=vcovHC(imp5.m1, type="HC0", cluster="group"))  
imp5.m1.coef
imp5.m1.frame=ctdf(imp5.m1.coef)

### Family & children
imp6.m1 <- plm(imp6.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                 c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", data)
imp6.m1.coef <- coeftest(imp6.m1, vcov=vcovHC(imp6.m1, type="HC0", cluster="group"))  
imp6.m1.coef
imp6.m1.frame=ctdf(imp6.m1.coef)

### Left-right ideology
lr.m1 <- plm(lr.s ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
               c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", data)
lr.m1.coef <- coeftest(lr.m1, vcov=vcovHC(lr.m1, type="HC0", cluster="group"))  
lr.m1.coef
lr.m1.frame=ctdf(lr.m1.coef)

################################################################################

### Figure S5: Effects of candidate disability on citizen perceptions of candidates’ 
# issue importance

t1.labels <- c("Belief", "Coefficient", "SE")
t1.row1 <- c("Social security & welfare", imp1.m1.frame[2,1], imp1.m1.frame[2,2] )
t1.row2 <- c("Military & defence", imp2.m1.frame[2,1], imp2.m1.frame[2,2] )
t1.row3 <- c("Healthcare", imp3.m1.frame[2,1], imp3.m1.frame[2,2] )
t1.row4 <- c("Minority rights", imp4.m1.frame[2,1], imp4.m1.frame[2,2] )
t1.row5 <- c("Economy", imp5.m1.frame[2,1], imp5.m1.frame[2,2] )
t1.row6 <- c("Family & children", imp6.m1.frame[2,1], imp6.m1.frame[2,2] )
t1.row7 <- c("Left-right", lr.m1.frame[2,1], lr.m1.frame[2,2] )

t1 <- data.frame(rbind(t1.row1,t1.row2,t1.row3,t1.row4,t1.row5,t1.row6,t1.row7))
names(t1) <- t1.labels
t1$Coefficient <- as.numeric(as.character(t1$Coefficient))
t1$SE <- as.numeric(as.character(t1$SE))
t1$Belief <- factor(t1$Belief, levels = c("Left-right", "Family & children", "Economy", 
                                          "Minority rights","Healthcare","Military & defence","Social security & welfare"))
interval <- 1.96  # 95% multiplier

plot <- ggplot(t1)
plot <- plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
plot <- plot + geom_pointrange(aes(x = Belief, y = Coefficient, ymin = Coefficient - SE*interval,
                                   ymax = Coefficient + SE*interval),
                               fill = "WHITE")
plot <- plot + ggtitle(" ") + theme_bw() +
  xlab(" ") +  ylab("Effect of candidate disability on perceived priorties and ideology") + 
  ylim(-0.08, 0.08) + coord_flip() 
print(plot)
ggsave("FigS5.png", dpi=300)

################################################################################

#### Table S3: Linear regression of perceived candidate preferences on candidate disability
  # (split into 2 tables as code doesn't run otherwise)

stargazer(imp1.m1.coef, imp2.m1.coef, imp3.m1.coef, imp4.m1.coef, 
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS3a.doc")

stargazer(imp5.m1.coef,
          imp6.m1.coef, lr.m1.coef,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS3b.doc")

# to get the N for each model: 
stargazer(imp1.m1, imp2.m1, imp3.m1, imp4.m1, imp5.m1,
          imp6.m1, lr.m1,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS3-n.doc")

################################################################################
################################################################################
################################################################################

#### Table S4 & Figure S6: Linear regression of citizen preferences on disability status and identity 

# Since these analyses are on the respondent level, we only use one observation per respondent, 
# by subsetting with candidate==1

data$disid3 <- relevel(as.factor(data$disid3), ref = "Non-disabled")

imp1own <- lm(imp.own1.n ~ disid3 + age + gender, data=subset(data, candidate==1))
frame1 <- summary(imp1own)
imp2own <- lm(imp.own2.n ~ disid3 + age + gender, data=subset(data, candidate==1))
frame2 <- summary(imp2own)
imp3own <- lm(imp.own3.n ~ disid3 + age + gender, data=subset(data, candidate==1))
frame3 <- summary(imp3own)
imp4own <- lm(imp.own4.n ~ disid3 + age + gender, data=subset(data, candidate==1))
frame4 <- summary(imp4own)
imp5own <- lm(imp.own5.n ~ disid3 + age + gender, data=subset(data, candidate==1))
frame5 <- summary(imp5own)
imp6own <- lm(imp.own6.n ~ disid3 + age + gender, data=subset(data, candidate==1))
frame6 <- summary(imp6own)
lrown <- lm(lrown.s ~ disid3 + age + gender, data=subset(data, candidate==1))
frame7 <- summary(lrown)

################################################################################

### Figure S6: Citizens’ issue importance (based on linear regressions, controlling for age and gender) 

t1.labels <- c("Belief", "Disability", "Coefficient", "SE")
t1.row1 <- c("Social security & welfare", "ID", frame1$coefficients[2,1], frame1$coefficients[2,2])
t1.row2 <- c("Military & defence", "ID", frame2$coefficients[2,1], frame2$coefficients[2,2])
t1.row3 <- c("Healthcare", "ID", frame3$coefficients[2,1], frame3$coefficients[2,2])
t1.row4 <- c("Minority rights", "ID", frame4$coefficients[2,1], frame4$coefficients[2,2])
t1.row5 <- c("Economy", "ID", frame5$coefficients[2,1], frame5$coefficients[2,2])
t1.row6 <- c("Family & children", "ID", frame6$coefficients[2,1], frame6$coefficients[2,2])
t1.row7 <- c("Left-right", "ID", frame7$coefficients[2,1], frame7$coefficients[2,2])
t1.row8 <- c("Social security & welfare", "No ID", frame1$coefficients[3,1], frame1$coefficients[3,2])
t1.row9 <- c("Military & defence", "No ID", frame2$coefficients[3,1], frame2$coefficients[3,2])
t1.row10 <- c("Healthcare", "No ID", frame3$coefficients[3,1], frame3$coefficients[3,2])
t1.row11 <- c("Minority rights", "No ID", frame4$coefficients[3,1], frame4$coefficients[3,2])
t1.row12 <- c("Economy", "No ID", frame5$coefficients[3,1], frame5$coefficients[3,2])
t1.row13 <- c("Family & children", "No ID", frame6$coefficients[3,1], frame6$coefficients[3,2])
t1.row14 <- c("Left-right", "No ID", frame7$coefficients[3,1], frame7$coefficients[3,2])

t1 <- data.frame(rbind(t1.row1,t1.row2,t1.row3,t1.row4,t1.row5,t1.row6,t1.row7,
                       t1.row8,t1.row9,t1.row10,t1.row11,t1.row12,t1.row13,t1.row14))
names(t1) <- t1.labels
t1$Coefficient <- as.numeric(as.character(t1$Coefficient))
t1$SE <- as.numeric(as.character(t1$SE))
t1$Belief <- factor(t1$Belief, levels = c("Left-right", "Family & children", "Economy", 
                                          "Minority rights","Healthcare","Military & defence","Social security & welfare"))
interval <- 1.96  # 95% multiplier

plot <- ggplot(t1)
plot <- plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
plot <- plot + geom_pointrange(aes(x = Belief, y = Coefficient, ymin = Coefficient - SE*interval,
                                   ymax = Coefficient + SE*interval, colour=Disability, shape=Disability),
                               lwd = 1/2, position = position_dodge(width = 1/2),
                               fill = "WHITE")
plot <- plot + ggtitle(" ") + theme_bw() +
  xlab(" ") +  ylab("Effect of citizen disability on priorities and ideology") + ylim(-0.08, 0.08) + coord_flip() + 
  scale_colour_manual(name="Citizen disability \n(reference category= \nnon-disabled)", values = c("black", "black")) +
  scale_shape_manual(guide="none", values=c(15,16)) +
  guides(colour = guide_legend(override.aes = list(shape = c(16,15)), reverse = TRUE))
print(plot)
ggsave("FigS6.png", dpi=300)

################################################################################

### Table S4: Linear regression of citizen preferences on disability status and identity

stargazer(imp1own, imp2own, imp3own, imp4own, imp5own, imp6own, lrown, 
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS4.doc")

################################################################################
################################################################################
################################################################################

##### SI 4: Mediation analysis 

# Make sure you've set the seed at the start to reproduce the results

### Table S5: Mediation models among disabled citizens with a disability identity

base.model <- lm(rep1.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                   c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country, 
                 data=subset(data, disid3=="ID" & !is.na(congboth)))
summary(base.model) 
med.fit <- lm(congboth ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country, 
              data=subset(data, disid3=="ID" & !is.na(rep1.nor)))
summary(med.fit)
out.fit.lm <- lm(rep1.nor ~ congboth + c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                   c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country,  
                 data=subset(data, disid3=="ID"))
summary(out.fit.lm)

med.out <- mediate(med.fit, out.fit.lm, treat = "c.disdummy", mediator = "congboth",
                   robustSE = T, sims = 1000)
summary(med.out)

# Sensitivity analysis (Figure S7: left)
sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out)
png(file="FigS7a.png")
plot(sens.out, sens.par = "rho", main = "Disability ID", ylim = c(-0.2, 0.2))
dev.off()

stargazer(base.model, med.fit, out.fit.lm,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS5.doc")

################################################################################

### Table S6: Mediation models among disabled citizens without a disability identity

base.model <- lm(rep1.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                   c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country, 
                 data=subset(data, disid3=="No ID" & !is.na(congboth)))
summary(base.model) 

med.fit <- lm(congboth ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country, 
              data=subset(data, disid3=="No ID" & !is.na(rep1.nor)))
summary(med.fit)
out.fit.lm <- lm(rep1.nor ~ congboth + c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                   c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country,  
                 data=subset(data, disid3=="No ID"))
summary(out.fit.lm)

med.out <- mediate(med.fit, out.fit.lm, treat = "c.disdummy", mediator = "congboth",
                   robustSE = T, sims = 1000)
summary(med.out)

# sensitivity analysis (Figure S7: middle)
sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out)
png(file="FigS7b.png")
plot(sens.out, sens.par = "rho", main = "No disability ID", ylim = c(-0.2, 0.2))
dev.off()

stargazer(base.model, med.fit, out.fit.lm,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="FigureS6.doc")

################################################################################

### Table S7: Mediation models among non-disabled citizens 

base.model <- lm(rep1.nor ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                   c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country, 
                 data=subset(data, disid3=="Non-disabled" & !is.na(congboth)))
summary(base.model) 

med.fit <- lm(congboth ~ c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country, 
              data=subset(data, disid3=="Non-disabled" & !is.na(rep1.nor)))
summary(med.fit)
out.fit.lm <- lm(rep1.nor ~ congboth + c.disdummy + c.age + c.female + c.minority + c.job + c.exp + 
                   c.child.num + c.heldoffice + c.party3 + as.factor(candidate) + country,  
                 data=subset(data, disid3=="Non-disabled"))
summary(out.fit.lm)

med.out <- mediate(med.fit, out.fit.lm, treat = "c.disdummy", mediator = "congboth",
                   robustSE = T, sims = 1000)
summary(med.out)

# Sensitivity analysis(Figure S7: right)
sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out)
png(file="FigS7c.png")
plot(sens.out, sens.par = "rho", main = "Non-disabled", ylim = c(-0.2, 0.2))
dev.off()

stargazer(base.model, med.fit, out.fit.lm, 
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS7.doc")

################################################################################

### Table S8: enter results from mediation analysis manually 

################################################################################
################################################################################
################################################################################

##### SI 5: Analysis of relationships between candidates’ and respondents’ disability 
# status and identity

# Since data is on the candidate level, we need to generate a variable that indicates 
# how many disabled candidates a respondent has seen: 

# dis_n: proportion of candidates seen who are disabled by respondent ID

data <- data %>%
  group_by(ID) %>%
  mutate(dis_n = mean(c.disdummy)) %>% ungroup
table(data$dis_n, useNA = "a")

# dis_b: binary variable indicating whether respondent saw disabled candidate

data$dis_b <- NA
data$dis_b[data$dis_n==0] <- 0 # no
data$dis_b[data$dis_n>0] <- 1 # yes
table(data$dis_b, useNA = "a")

# For the analysis, we use only one observation per respondent, therefore we create a
# data subset that only includes rows for candidate=1:
subdata <- subset(data, candidate==1)

################################################################################

### Table S9: Relationship between candidate disability and respondent disability status 

CrossTable(y=subdata$disown, x=subdata$dis_b, expected=T, prop.r=T, prop.c=F, prop.t=F,
           prop.chisq=F, chisq=T, format="SPSS")

# t-test testing whether the mean number of disabled candidates seen varies with 
# respondents' own disability status (reported in text)
t.test(dis_n ~ disown, data=subdata)

################################################################################

### Table S10: Relationship between candidate disability and respondent disability group identity 

CrossTable(y=subdata$disid, x=subdata$dis_b, expected=T, prop.r=T, prop.c=F, prop.t=F,
           prop.chisq=F, chisq=T, format="SPSS")

# t-test testing whether the mean number of disabled candidates seen varies with 
# respondents' own disability identity (reported in text)
t.test(dis_n ~ disid, data=subdata)

################################################################################

### Table S11: Linear regression of representation perceptions on candidate disability 
# (Model 2, Table S2) among respondents who saw at least one disabled candidate 

data$disid3 <- as.character(data$disid3)
data$disid3 <- as.factor(data$disid3)

rep1.m2.saw <- plm(rep1.nor ~ c.disdummy*relevel(disid3, ref="Non-disabled") + c.age + c.female + c.minority + c.job + c.exp + 
                     c.child.num + c.heldoffice + c.party3 + country + as.factor(candidate), model="p", index="ID", 
                   data=subset(data, dis_b==1))
rep1.m2.saw.coef <- coeftest(rep1.m2.saw, vcov=vcovHC(rep1.m2.saw, type="HC0", cluster="group"))  
rep1.m2.saw.coef

stargazer(rep1.m2.saw,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="html", out="TableS11.doc")



